knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
report_hurdle = TRUE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE # ! set to TRUE !

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1635 0.0000 5.0000 241

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Own Actionplan",
    'Partner Actionplan',
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,10),
  "Between-Person Effects" = c(11,17),
  "Random Effects" = c(18, 24), 
  "Additional Parameters" = c(25,25)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,10+5),
  "Between-Person Effects" = c(11+5,17+5),
  "Random Effects" = c(18+5, 24+5), 
  "Additional Parameters" = c(25+5,25+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_plan_selfPlan',
    'hu_plan_partnerPlan',
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_fixed_hu_ordinal <- c(
  model_rows_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed_hu[3:length(model_rows_fixed_hu)]
)


model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_hu_ordinal <- c(model_rows_random_hu,'disc')
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Own actionplan",
    'Partner actionplan', 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Own actionplan",
    'Hu Partner actionplan', 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )



model_rownames_fixed_hu_ordinal <- c(
  model_rownames_fixed_hu[1:2],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed_hu[3:length(model_rownames_fixed_hu)]
)



model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_hu_ordinal <- c(model_rownames_random_hu,'disc')
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,12),
  "Conditional Between-Person Effects" = c(13,19),
  
  "Hurdle Within-Person Effects" = c(20,29),
  "Hurdle Between-Person Effects" = c(30,36),
  
  "Random Effects" = c(37, 50), 
  "Additional Parameters" = c(51,51)
  )

rows_to_pack_hu_ordinal <- list(
  "Conditional Within-Person Effects" = c(3+5,12+5),
  "Conditional Between-Person Effects" = c(13+5,19+5),
  
  "Hurdle Within-Person Effects" = c(20+5,29+5),
  "Hurdle Between-Person Effects" = c(30+5,36+5),
  
  "Random Effects" = c(37+5, 50+5), 
  "Additional Parameters" = c(51+5,51+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 7, observations = 3736, p-value = 0.56
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000267666 0.002817184
## sample estimates:
## outlier frequency (expected: 0.00149892933618844 ) 
##                                        0.001873662
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# summarize with rope range for hurdle part
summary_pa_sub_hurdle <- summarize_brms(
  pa_sub, 
  stats_to_report = c('pd','ROPE'),
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Warning in summarize_brms(pa_sub, stats_to_report = c("pd", "ROPE"), rope_range
## = c(-0.18, : Coefficients were exponentiated. Double check if this was
## intended.
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub_continuous <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report, rope_range
## = rope_range_continuous, : Coefficients were exponentiated. Double check if
## this was intended.
# Replace only the ROPE and % in Rope columns for rows with 'Hu'
summary_pa_sub <- summary_pa_sub_continuous

columns_to_replace <- c("ROPE", "inside ROPE")

summary_pa_sub[grepl('Hu', rownames(summary_pa_sub)), columns_to_replace] <- 
  summary_pa_sub_hurdle[grepl('Hu', rownames(summary_pa_sub_hurdle)), columns_to_replace]

# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack_hu)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 37.92*** 2.64 [33.02, 43.54] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 12653 22304
Hurdle Intercept 0.26*** 0.04 [ 0.19, 0.35] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 10801 19932
Conditional Within-Person Effects
Daily persuasion experienced 1.03 0.03 [ 0.98, 1.09] 0.859 [0.92, 1.08] 0.967 0.019 Very Strong Evidence for Null 1.000 17362 24808
Daily persuasion utilized (partner’s view) 1.03 0.02 [ 0.99, 1.08] 0.916 [0.92, 1.08] 0.978 0.023 Very Strong Evidence for Null 1.000 22959 28221
Daily pressure experienced 0.91* 0.04 [ 0.82, 1.00] 0.975 [0.92, 1.08] 0.343 0.154 Moderate Evidence for Null 1.000 36639 29172
Daily pressure utilized (partner’s view) 0.95 0.04 [ 0.86, 1.03] 0.894 [0.92, 1.08] 0.700 0.038 Strong Evidence for Null 1.000 37501 28659
Daily pushing experienced 0.99 0.03 [ 0.92, 1.05] 0.674 [0.92, 1.08] 0.966 0.014 Very Strong Evidence for Null 1.000 31510 30242
Daily pushing utilized (partner’s view) 0.96 0.03 [ 0.90, 1.02] 0.914 [0.92, 1.08] 0.882 0.031 Strong Evidence for Null 1.000 33140 30422
Day 0.99 0.06 [ 0.88, 1.11] 0.577 [0.92, 1.08] 0.790 0.025 Very Strong Evidence for Null 1.000 50335 29874
Own actionplan 1.33*** 0.06 [ 1.21, 1.45] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 46130 29477
Partner actionplan 1.08 0.05 [ 0.99, 1.17] 0.962 [0.92, 1.08] 0.518 0.084 Strong Evidence for Null 1.000 47984 32207
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.03 0.16 [ 0.77, 1.40] 0.590 [0.92, 1.08] 0.386 0.061 Strong Evidence for Null 1.000 7638 14348
Mean persuasion utilized (partner’s view) 0.99 0.15 [ 0.73, 1.34] 0.524 [0.92, 1.08] 0.396 0.061 Strong Evidence for Null 1.000 7814 14204
Mean pressure experienced 1.18 0.21 [ 0.83, 1.67] 0.826 [0.92, 1.08] 0.227 0.110 Moderate Evidence for Null 1.000 11628 20761
Mean pressure utilized (partner’s view) 0.92 0.17 [ 0.65, 1.32] 0.671 [0.92, 1.08] 0.300 0.081 Strong Evidence for Null 1.000 12201 21614
Mean pushing experienced 1.21 0.27 [ 0.77, 1.88] 0.797 [0.92, 1.08] 0.193 0.128 Moderate Evidence for Null 1.000 11020 18100
Mean pushing utilized (partner’s view) 1.30 0.30 [ 0.82, 2.06] 0.871 [0.92, 1.08] 0.141 0.175 Moderate Evidence for Null 1.000 11072 19999
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.59*** 0.11 [ 1.39, 1.85] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 26190 27004
Hu Daily persuasion utilized (partner’s view) 1.34*** 0.09 [ 1.19, 1.54] 1.000 [0.84, 1.20] 0.035 >100 Overwhelming Evidence 1.000 30919 27774
Hu Daily pressure experienced 0.95 0.14 [ 0.69, 1.29] 0.622 [0.84, 1.20] 0.742 0.079 Strong Evidence for Null 1.000 35225 27790
Hu Daily pressure utilized (partner’s view) 1.49* 0.28 [ 1.05, 2.29] 0.987 [0.84, 1.20] 0.113 1.098 Weak Evidence 1.000 34454 22964
Hu Daily pushing experienced 0.94 0.14 [ 0.71, 1.28] 0.665 [0.84, 1.20] 0.735 0.083 Strong Evidence for Null 1.000 23390 25680
Hu Daily pushing utilized (partner’s view) 1.28* 0.15 [ 1.04, 1.64] 0.990 [0.84, 1.20] 0.264 0.819 Weak Evidence for Null 1.000 35469 27582
Hu Day 0.86 0.13 [ 0.64, 1.15] 0.843 [0.84, 1.20] 0.566 0.124 Moderate Evidence for Null 1.000 61855 30632
Hu Own actionplan 9.48*** 0.97 [ 7.77, 11.60] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 45197 31239
Hu Partner actionplan 1.17 0.12 [ 0.96, 1.42] 0.941 [0.84, 1.20] 0.583 0.178 Moderate Evidence for Null 1.000 44112 30129
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.40 0.50 [ 0.69, 2.84] 0.825 [0.84, 1.20] 0.256 0.280 Moderate Evidence for Null 1.001 8405 15769
Hu Mean persuasion utilized (partner’s view) 1.33 0.47 [ 0.65, 2.67] 0.783 [0.84, 1.20] 0.289 0.240 Moderate Evidence for Null 1.000 8466 14879
Hu Mean pressure experienced 0.40* 0.17 [ 0.18, 0.90] 0.986 [0.84, 1.20] 0.034 2.480 Weak Evidence 1.000 14450 24340
Hu Mean pressure utilized (partner’s view) 0.51 0.21 [ 0.22, 1.15] 0.948 [0.84, 1.20] 0.096 0.776 Weak Evidence for Null 1.000 14670 22245
Hu Mean pushing experienced 1.17 0.62 [ 0.41, 3.28] 0.616 [0.84, 1.20] 0.252 0.278 Moderate Evidence for Null 1.000 12680 20855
Hu Mean pushing utilized (partner’s view) 2.07 1.10 [ 0.71, 5.84] 0.916 [0.84, 1.20] 0.105 0.687 Weak Evidence for Null 1.000 12400 20808
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 0.04 [0.23, 0.41] NA NA NA NA NA 1.000 10982 19107
sd(Hurdle Intercept) 0.73 0.10 [0.56, 0.99] NA NA NA NA NA 1.000 10861 19564
sd(Daily persuasion experienced) 0.11 0.02 [0.07, 0.17] NA NA NA NA NA 1.000 20621 28034
sd(Daily persuasion utilized (partner’s view)) 0.08 0.02 [0.04, 0.13] NA NA NA NA NA 1.000 21526 23726
sd(Daily pressure experienced) 0.06 0.05 [0.00, 0.22] NA NA NA NA NA 1.000 16670 18381
sd(Daily pressure utilized (partner’s view)) 0.06 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 17699 15819
sd(Daily pushing experienced) 0.09 0.04 [0.01, 0.17] NA NA NA NA NA 1.000 11617 9522
sd(Daily pushing utilized (partner’s view)) 0.08 0.03 [0.01, 0.15] NA NA NA NA NA 1.000 12829 10346
sd(Hu Daily persuasion experienced) 0.21 0.09 [0.03, 0.40] NA NA NA NA NA 1.001 10174 8865
sd(Hu Daily persuasion utilized (partner’s view)) 0.19 0.09 [0.02, 0.37] NA NA NA NA NA 1.000 9864 8686
sd(Hu Daily pressure experienced) 0.18 0.16 [0.01, 0.68] NA NA NA NA NA 1.000 16735 17846
sd(Hu Daily pressure utilized (partner’s view)) 0.23 0.21 [0.01, 0.89] NA NA NA NA NA 1.000 16584 19362
sd(Hu Daily pushing experienced) 0.57 0.16 [0.30, 0.98] NA NA NA NA NA 1.000 18443 25123
sd(Hu Daily pushing utilized (partner’s view)) 0.22 0.14 [0.02, 0.54] NA NA NA NA NA 1.000 13769 12571
Additional Parameters
sigma 0.68 0.01 [0.65, 0.70] NA NA NA NA NA 1.000 54627 28168 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.76), b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.74). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## This is posterior version 1.6.0
## 
## Attaching package: 'posterior'
## The following object is masked from 'package:bayesplot':
## 
##     rhat
## The following objects are masked from 'package:stats':
## 
##     mad, sd, var
## The following objects are masked from 'package:base':
## 
##     %in%, match
## Registering fonts with R
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 80 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00828

$persuasion_partner_cw

## Warning: Removed 61 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00653

$pressure_self_cw

## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.011

$pressure_partner_cw

## Picking joint bandwidth of 0.0184

$pushing_self_cw

## Picking joint bandwidth of 0.0103

$pushing_partner_cw

## Picking joint bandwidth of 0.0101

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  )

for (i in 1:length(conds_eff)) {
  effname <- names(conds_eff)[i]
  eff_plot <- conds_eff[[i]]
  x_label_i <- x_label[[i]]
  rmarkdown::render(
    "C:/Users/kueng/DataAnalysis/02TimeAndTiesControl/Output/Plots/BeautifulPlotWithNote.Rmd", 
    output_file = paste0('C:/Users/kueng/DataAnalysis/02TimeAndTiesControl/Output/Plots/Graphic_', effname, '.pdf'),
    params = list(
      p_i = eff_plot,
      p_name = effname,
      x_label = x_label_i
      ),
    envir = new.env(),
    quiet = TRUE
  )
}

print('done')

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 28, observations = 3337, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.001198681 0.004794726
## sample estimates:
## outlier frequency (expected: 0.00297872340425532 ) 
##                                         0.00839077
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 111.17*** 5.93 [99.95, 123.93] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.001 6903 13661
Within-Person Effects
Daily persuasion experienced 1.03 0.02 [ 1.00, 1.06] 0.953 [0.94, 1.07] 0.995 0.026 Very Strong Evidence for Null 1.000 26993 29364
Daily persuasion utilized (partner’s view) 1.02 0.02 [ 0.99, 1.05] 0.864 [0.94, 1.07] 0.998 0.011 Very Strong Evidence for Null 1.000 31647 30752
Daily pressure experienced 0.95 0.03 [ 0.88, 1.01] 0.947 [0.94, 1.07] 0.602 0.055 Strong Evidence for Null 1.000 47423 30997
Daily pressure utilized (partner’s view) 0.98 0.03 [ 0.92, 1.05] 0.694 [0.94, 1.07] 0.916 0.015 Very Strong Evidence for Null 1.000 50157 31658
Daily pushing experienced 1.02 0.02 [ 0.96, 1.07] 0.731 [0.94, 1.07] 0.974 0.012 Very Strong Evidence for Null 1.000 35878 29732
Daily pushing utilized (partner’s view) 1.00 0.02 [ 0.96, 1.04] 0.510 [0.94, 1.07] 0.997 0.008 Very Strong Evidence for Null 1.000 48187 32131
Day 0.97 0.03 [ 0.91, 1.04] 0.788 [0.94, 1.07] 0.845 0.020 Very Strong Evidence for Null 1.000 77038 29175
Own Actionplan 1.06* 0.03 [ 1.01, 1.12] 0.993 [0.94, 1.07] 0.530 0.225 Moderate Evidence for Null 1.000 57658 31464
Partner Actionplan 1.05 0.03 [ 1.00, 1.10] 0.972 [0.94, 1.07] 0.753 0.060 Strong Evidence for Null 1.000 57893 31743
Between-Person Effects
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 0.053 Strong Evidence for Null 1.000 41154 26326
Mean persuasion experienced 1.10 0.16 [ 0.83, 1.46] 0.756 [0.94, 1.07] 0.279 0.072 Strong Evidence for Null 1.001 5287 11364
Mean persuasion utilized (partner’s view) 0.98 0.14 [ 0.73, 1.30] 0.562 [0.94, 1.07] 0.344 0.059 Strong Evidence for Null 1.001 5210 10807
Mean pressure experienced 0.99 0.14 [ 0.74, 1.33] 0.531 [0.94, 1.07] 0.344 0.058 Strong Evidence for Null 1.001 7594 15835
Mean pressure utilized (partner’s view) 0.97 0.14 [ 0.74, 1.29] 0.577 [0.94, 1.07] 0.347 0.058 Strong Evidence for Null 1.001 6847 13727
Mean pushing experienced 0.95 0.19 [ 0.63, 1.43] 0.600 [0.94, 1.07] 0.241 0.084 Strong Evidence for Null 1.001 8398 15531
Mean pushing utilized (partner’s view) 1.22 0.25 [ 0.82, 1.84] 0.840 [0.94, 1.07] 0.156 0.133 Moderate Evidence for Null 1.001 8149 15791
Random Effects
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.916 [0.94, 1.07] 1.000 0.000 Very Strong Evidence for Null 1.000 50144 33445
sd(Intercept) 0.29 0.04 [0.23, 0.39] NA NA NA NA NA 1.001 8384 17687
sd(Daily persuasion experienced) 0.05 0.01 [0.02, 0.08] NA NA NA NA NA 1.000 16704 10988
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 [0.03, 0.09] NA NA NA NA NA 1.000 19719 20938
sd(Daily pressure experienced) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 19449 21382
sd(Daily pressure utilized (partner’s view)) 0.03 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 23291 19151
sd(Daily pushing experienced) 0.07 0.04 [0.01, 0.15] NA NA NA NA NA 1.000 9926 12203
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.03 0.03 [0.00, 0.09] NA NA NA NA NA 1.000 16153 18099
sigma 0.57 0.01 [0.56, 0.59] NA NA NA NA NA 1.000 65401 28026 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.89), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.74), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.72), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.77), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 30, observations = 3736, p-value =
## 4.114e-10
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.005424162 0.011443600
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.008029979
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.67*** 0.10 [ 3.46, 3.88] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 4513 9768
Within-Person Effects
Daily persuasion experienced 0.00 0.02 [-0.04, 0.04] 0.523 [-0.11, 0.11] 1.000 0.004 Very Strong Evidence for Null 1.000 35272 26862
Daily persuasion utilized (partner’s view) 0.02 0.02 [-0.02, 0.07] 0.810 [-0.11, 0.11] 1.000 0.007 Very Strong Evidence for Null 1.000 30965 30203
Daily pressure experienced -0.03 0.05 [-0.14, 0.07] 0.717 [-0.11, 0.11] 0.937 0.012 Very Strong Evidence for Null 1.000 36650 27776
Daily pressure utilized (partner’s view) -0.03 0.05 [-0.15, 0.08] 0.694 [-0.11, 0.11] 0.930 0.012 Very Strong Evidence for Null 1.000 34035 26217
Daily pushing experienced 0.01 0.03 [-0.06, 0.07] 0.577 [-0.11, 0.11] 0.999 0.007 Very Strong Evidence for Null 1.000 42746 31073
Daily pushing utilized (partner’s view) 0.07* 0.03 [ 0.00, 0.14] 0.976 [-0.11, 0.11] 0.895 0.058 Strong Evidence for Null 1.000 31862 28826
Day 0.26*** 0.06 [ 0.15, 0.37] 1.000 [-0.11, 0.11] 0.005 >100 Overwhelming Evidence 1.000 58650 29078
Own Actionplan 0.10** 0.04 [ 0.03, 0.18] 0.996 [-0.11, 0.11] 0.618 0.268 Moderate Evidence for Null 1.000 48149 31839
Partner Actionplan -0.03 0.04 [-0.11, 0.05] 0.792 [-0.11, 0.11] 0.982 0.011 Very Strong Evidence for Null 1.000 47941 31849
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Mean persuasion experienced 0.34 0.27 [-0.20, 0.90] 0.892 [-0.11, 0.11] 0.154 0.120 Moderate Evidence for Null 1.001 3933 7303
Mean persuasion utilized (partner’s view) 0.22 0.27 [-0.32, 0.78] 0.794 [-0.11, 0.11] 0.236 0.076 Strong Evidence for Null 1.001 3906 7095
Mean pressure experienced -0.30 0.27 [-0.85, 0.24] 0.866 [-0.11, 0.11] 0.181 0.102 Moderate Evidence for Null 1.001 4698 9345
Mean pressure utilized (partner’s view) -0.31 0.27 [-0.86, 0.22] 0.877 [-0.11, 0.11] 0.170 0.110 Moderate Evidence for Null 1.001 4686 9233
Mean pushing experienced 0.21 0.38 [-0.55, 0.97] 0.707 [-0.11, 0.11] 0.204 0.087 Strong Evidence for Null 1.000 6174 12395
Mean pushing utilized (partner’s view) 0.37 0.38 [-0.39, 1.13] 0.835 [-0.11, 0.11] 0.149 0.124 Moderate Evidence for Null 1.000 6234 12509
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
sd(Intercept) 0.60 0.07 [0.48, 0.78] NA NA NA NA NA 1.000 7398 14380
sd(Daily persuasion experienced) 0.04 0.03 [0.00, 0.10] NA NA NA NA NA 1.000 10459 14097
sd(Daily persuasion utilized (partner’s view)) 0.07 0.03 [0.01, 0.13] NA NA NA NA NA 1.001 11847 10659
sd(Daily pressure experienced) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 16312 19973
sd(Daily pressure utilized (partner’s view)) 0.08 0.07 [0.00, 0.26] NA NA NA NA NA 1.000 15919 18986
sd(Daily pushing experienced) 0.05 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 15239 17528
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.07 0.05 [0.00, 0.17] NA NA NA NA NA 1.000 13634 14941
sigma 0.96 0.01 [0.94, 0.98] NA NA NA NA NA 1.000 54740 29154 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.82), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.8), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.8), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.83), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.77), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.89). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 4 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.001322751
## sample estimates:
## outlier frequency (expected: 0.000317460317460317 ) 
##                                         0.002645503
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercepts
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.36*** 1.02 [ 1.88, 6.23] 1.000 [0.84, 1.20] 0.000 69.572 Very Strong Evidence 1.000 33161 30736
Intercept[2] 7.31*** 2.30 [ 4.00, 13.82] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 33914 29999
Intercept[3] 20.38*** 6.87 [ 10.74, 40.28] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 35591 31079
Intercept[4] 89.21*** 34.75 [ 42.90, 195.22] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 39807 32159
Intercept[5] 3029.60*** 1996.62 [916.76, 12308.87] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 49408 30767
Within-Person Effects
Daily persuasion experienced 0.84* 0.07 [ 0.71, 0.99] 0.982 [0.84, 1.20] 0.549 0.297 Moderate Evidence for Null 1.000 33919 28293
Daily persuasion utilized (partner’s view) 1.02 0.10 [ 0.83, 1.23] 0.583 [0.84, 1.20] 0.925 0.040 Strong Evidence for Null 1.000 32260 28470
Daily pressure experienced 1.84* 0.37 [ 1.16, 2.66] 0.992 [0.84, 1.20] 0.031 2.539 Weak Evidence 1.000 21582 22822
Daily pressure utilized (partner’s view) 1.24 0.30 [ 0.69, 2.11] 0.807 [0.84, 1.20] 0.376 0.158 Moderate Evidence for Null 1.000 23305 19657
Daily pushing experienced 1.22 0.13 [ 0.98, 1.52] 0.962 [0.84, 1.20] 0.445 0.217 Moderate Evidence for Null 1.000 28606 27455
Daily pushing utilized (partner’s view) 0.94 0.12 [ 0.71, 1.22] 0.693 [0.84, 1.20] 0.776 0.059 Strong Evidence for Null 1.000 32445 27556
Day 1.48 0.51 [ 0.75, 2.88] 0.872 [0.84, 1.20] 0.221 0.266 Moderate Evidence for Null 1.000 47791 31765
Own Actionplan 0.85 0.24 [ 0.49, 1.49] 0.714 [0.84, 1.20] 0.410 0.134 Moderate Evidence for Null 1.000 35305 32306
Partner Actionplan 0.92 0.24 [ 0.55, 1.52] 0.632 [0.84, 1.20] 0.493 0.110 Moderate Evidence for Null 1.000 37423 31789
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Mean persuasion experienced 1.10 0.56 [ 0.39, 3.02] 0.572 [0.84, 1.20] 0.269 0.211 Moderate Evidence for Null 1.000 15958 22898
Mean persuasion utilized (partner’s view) 1.37 0.77 [ 0.45, 4.27] 0.709 [0.84, 1.20] 0.219 0.265 Moderate Evidence for Null 1.000 16443 22952
Mean pressure experienced 3.52* 1.93 [ 1.22, 10.68] 0.990 [0.84, 1.20] 0.020 3.211 Moderate Evidence 1.000 18630 25099
Mean pressure utilized (partner’s view) 1.16 0.67 [ 0.36, 3.54] 0.602 [0.84, 1.20] 0.235 0.236 Moderate Evidence for Null 1.000 17309 24737
Mean pushing experienced 1.25 0.93 [ 0.30, 5.72] 0.623 [0.84, 1.20] 0.187 0.310 Weak Evidence for Null 1.000 20161 26058
Mean pushing utilized (partner’s view) 0.11* 0.10 [ 0.02, 0.65] 0.993 [0.84, 1.20] 0.009 7.349 Moderate Evidence 1.000 26207 29097
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
sd(Intercept) 0.80 0.20 [0.46, 1.27] NA NA NA NA NA 1.000 15427 22186
sd(Daily persuasion experienced) 0.16 0.12 [0.01, 0.43] NA NA NA NA NA 1.000 7784 15596
sd(Daily persuasion utilized (partner’s view)) 0.21 0.14 [0.01, 0.52] NA NA NA NA NA 1.001 10675 13886
sd(Daily pressure experienced) 0.56 0.25 [0.11, 1.18] NA NA NA NA NA 1.000 8950 7907
sd(Daily pressure utilized (partner’s view)) 0.42 0.40 [0.02, 1.56] NA NA NA NA NA 1.001 9732 18185
sd(Daily pushing experienced) 0.21 0.14 [0.01, 0.51] NA NA NA NA NA 1.000 11006 14481
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.15 0.14 [0.01, 0.62] NA NA NA NA NA 1.000 17161 19837
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.8), b_Intercept[4] and b_Intercept[3] (r = 0.86), b_pressure_self_cb
##   and b_persuasion_self_cb (r = 0.71), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.79). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="01_FinalModelsPlan_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 35 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 2, observations = 756, p-value = 0.6663
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0003205434 0.0095234999
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002645503
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Sampling priors, please wait...
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.34*** 0.11 [0.17, 0.64] 1.000 [0.83, 1.20] 0.003 12.437 Strong Evidence 1.000 28741 31565
Within-Person Effects
Daily persuasion experienced 0.84 0.08 [0.68, 1.01] 0.968 [0.83, 1.20] 0.527 0.215 Moderate Evidence for Null 1.000 26792 22773
Daily persuasion utilized (partner’s view) 1.12 0.16 [0.84, 1.53] 0.780 [0.83, 1.20] 0.666 0.074 Strong Evidence for Null 1.000 20781 21793
Daily pressure experienced 1.97* 0.64 [1.01, 4.43] 0.976 [0.83, 1.20] 0.053 1.201 Weak Evidence 1.000 16708 16138
Daily pressure utilized (partner’s view) 1.42 0.60 [0.58, 4.18] 0.804 [0.83, 1.20] 0.234 0.241 Moderate Evidence for Null 1.000 16914 17573
Daily pushing experienced 1.33* 0.18 [1.03, 1.75] 0.986 [0.83, 1.20] 0.208 0.603 Weak Evidence for Null 1.000 24947 25401
Daily pushing utilized (partner’s view) 0.92 0.17 [0.62, 1.36] 0.677 [0.83, 1.20] 0.605 0.087 Strong Evidence for Null 1.000 28545 25800
Day 1.66 0.65 [0.77, 3.60] 0.905 [0.83, 1.20] 0.160 0.368 Weak Evidence for Null 1.000 39966 30260
Own Actionplan 0.87 0.28 [0.46, 1.64] 0.665 [0.83, 1.20] 0.395 0.139 Moderate Evidence for Null 1.000 30891 30274
Partner Actionplan 0.87 0.26 [0.49, 1.54] 0.681 [0.83, 1.20] 0.421 0.130 Moderate Evidence for Null 1.000 35248 29799
Between-Person Effects
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Mean persuasion experienced 1.90 1.16 [0.57, 6.60] 0.856 [0.83, 1.20] 0.135 0.430 Weak Evidence for Null 1.000 14375 21061
Mean persuasion utilized (partner’s view) 1.83 1.23 [0.50, 7.20] 0.819 [0.83, 1.20] 0.144 0.395 Weak Evidence for Null 1.000 14643 22514
Mean pressure experienced 18.34** 19.65 [2.51, 165.74] 0.998 [0.83, 1.20] 0.002 29.292 Strong Evidence 1.000 18604 23442
Mean pressure utilized (partner’s view) 2.29 2.54 [0.23, 19.25] 0.769 [0.83, 1.20] 0.095 0.613 Weak Evidence for Null 1.000 17099 24008
Mean pushing experienced 0.86 0.85 [0.12, 6.20] 0.562 [0.83, 1.20] 0.144 0.401 Weak Evidence for Null 1.000 16631 23514
Mean pushing utilized (partner’s view) 0.08* 0.09 [0.01, 0.68] 0.990 [0.83, 1.20] 0.009 6.280 Moderate Evidence 1.000 21332 27410
Random Effects
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
sd(Intercept) 1.15 0.25 [0.74, 1.76] NA NA NA NA NA 1.000 12310 20831
sd(Daily persuasion experienced) 0.22 0.14 [0.01, 0.51] NA NA NA NA NA 1.000 7209 12307
sd(Daily persuasion utilized (partner’s view)) 0.49 0.20 [0.13, 0.98] NA NA NA NA NA 1.000 11110 10788
sd(Daily pressure experienced) 1.08 0.55 [0.14, 2.43] NA NA NA NA NA 1.001 6642 6944
sd(Daily pressure utilized (partner’s view)) 0.83 0.68 [0.04, 2.75] NA NA NA NA NA 1.000 10859 14383
sd(Daily pushing experienced) 0.25 0.16 [0.02, 0.62] NA NA NA NA NA 1.000 10817 11169
Additional Parameters
sd(Daily pushing utilized (partner’s view)) 0.26 0.23 [0.01, 0.96] NA NA NA NA NA 1.000 13653 16584
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="01_FinalModelsPlan_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.41       0.4    -0.23     1.08       6.09
##   Post.Prob Star
## 1      0.86     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

process_model_component <- function(obj) {
  # Convert to string, modify, and evaluate
  name <- deparse(substitute(obj))
  if (report_hurdle) name <- paste0(name, '_hu')
  if (report_ordinal) name <- paste0(name, '_ordinal')
  return(get(name, envir = parent.frame()))
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  reactance_ordinal,
  is_reactance,
  
  stats_to_report = c('CI', 'pd'),
  
  model_rows_random = process_model_component(model_rows_random),
  model_rows_fixed = process_model_component(model_rows_fixed),
  model_rownames_random = process_model_component(model_rownames_random),
  model_rownames_fixed = process_model_component(model_rownames_fixed)
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “reactance_ordinal” [1] “is_reactance”

# pretty printing

summary_all_models <- all_models %>%
  print_df(rows_to_pack = process_model_component(rows_to_pack)) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = (length(all_models) / 5),  
      "Device-Based MVPA Log (Gaussian)" = (length(all_models) / 5), 
      "Mood Gaussian" = (length(all_models) / 5),
      "Reactance Ordinal" = (length(all_models) / 5),
      "Reactance Dichotome" = (length(all_models) / 5))
  )


export_xlsx(
  summary_all_models, 
  rows_to_pack = process_model_component(rows_to_pack),
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'header', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
print(summary_all_models)
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Ordinal
Reactance Dichotome
exp(Est.) pa_sub 95% CI pa_sub pd pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log pd pa_obj_log Est. mood_gauss 95% CI mood_gauss pd mood_gauss OR reactance_ordinal 95% CI reactance_ordinal pd reactance_ordinal OR is_reactance 95% CI is_reactance pd is_reactance
Intercept 37.92*** [33.02, 43.54] 1.000 111.17*** [99.95, 123.93] 1.000 3.67*** [ 3.46, 3.88] 1.000 NA NA NA 0.34*** [0.17, 0.64] 1.000
Hurdle Intercept 0.26*** [ 0.19, 0.35] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 1.03 [ 0.98, 1.09] 0.859 1.03 [ 1.00, 1.06] 0.953 0.00 [-0.04, 0.04] 0.523 0.84* [ 0.71, 0.99] 0.982 0.84 [0.68, 1.01] 0.968
Daily persuasion utilized (partner’s view) 1.03 [ 0.99, 1.08] 0.916 1.02 [ 0.99, 1.05] 0.864 0.02 [-0.02, 0.07] 0.810 1.02 [ 0.83, 1.23] 0.583 1.12 [0.84, 1.53] 0.780
Daily pressure experienced 0.91* [ 0.82, 1.00] 0.975 0.95 [ 0.88, 1.01] 0.947 -0.03 [-0.14, 0.07] 0.717 1.84* [ 1.16, 2.66] 0.992 1.97* [1.01, 4.43] 0.976
Daily pressure utilized (partner’s view) 0.95 [ 0.86, 1.03] 0.894 0.98 [ 0.92, 1.05] 0.694 -0.03 [-0.15, 0.08] 0.694 1.24 [ 0.69, 2.11] 0.807 1.42 [0.58, 4.18] 0.804
Daily pushing experienced 0.99 [ 0.92, 1.05] 0.674 1.02 [ 0.96, 1.07] 0.731 0.01 [-0.06, 0.07] 0.577 1.22 [ 0.98, 1.52] 0.962 1.33* [1.03, 1.75] 0.986
Daily pushing utilized (partner’s view) 0.96 [ 0.90, 1.02] 0.914 1.00 [ 0.96, 1.04] 0.510 0.07* [ 0.00, 0.14] 0.976 0.94 [ 0.71, 1.22] 0.693 0.92 [0.62, 1.36] 0.677
Day 0.99 [ 0.88, 1.11] 0.577 0.97 [ 0.91, 1.04] 0.788 0.26*** [ 0.15, 0.37] 1.000 1.48 [ 0.75, 2.88] 0.872 1.66 [0.77, 3.60] 0.905
Own actionplan 1.33*** [ 1.21, 1.45] 1.000 1.06* [ 1.01, 1.12] 0.993 0.10** [ 0.03, 0.18] 0.996 0.85 [ 0.49, 1.49] 0.714 0.87 [0.46, 1.64] 0.665
Partner actionplan 1.08 [ 0.99, 1.17] 0.962 1.05 [ 1.00, 1.10] 0.972 -0.03 [-0.11, 0.05] 0.792 0.92 [ 0.55, 1.52] 0.632 0.87 [0.49, 1.54] 0.681
Daily weartime NA NA NA 1.00*** [ 1.00, 1.00] 1.000 NA NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.03 [ 0.77, 1.40] 0.590 1.10 [ 0.83, 1.46] 0.756 0.34 [-0.20, 0.90] 0.892 1.10 [ 0.39, 3.02] 0.572 1.90 [0.57, 6.60] 0.856
Mean persuasion utilized (partner’s view) 0.99 [ 0.73, 1.34] 0.524 0.98 [ 0.73, 1.30] 0.562 0.22 [-0.32, 0.78] 0.794 1.37 [ 0.45, 4.27] 0.709 1.83 [0.50, 7.20] 0.819
Mean pressure experienced 1.18 [ 0.83, 1.67] 0.826 0.99 [ 0.74, 1.33] 0.531 -0.30 [-0.85, 0.24] 0.866 3.52* [ 1.22, 10.68] 0.990 18.34** [2.51, 165.74] 0.998
Mean pressure utilized (partner’s view) 0.92 [ 0.65, 1.32] 0.671 0.97 [ 0.74, 1.29] 0.577 -0.31 [-0.86, 0.22] 0.877 1.16 [ 0.36, 3.54] 0.602 2.29 [0.23, 19.25] 0.769
Mean pushing experienced 1.21 [ 0.77, 1.88] 0.797 0.95 [ 0.63, 1.43] 0.600 0.21 [-0.55, 0.97] 0.707 1.25 [ 0.30, 5.72] 0.623 0.86 [0.12, 6.20] 0.562
Mean pushing utilized (partner’s view) 1.30 [ 0.82, 2.06] 0.871 1.22 [ 0.82, 1.84] 0.840 0.37 [-0.39, 1.13] 0.835 0.11* [ 0.02, 0.65] 0.993 0.08* [0.01, 0.68] 0.990
Mean weartime NA NA NA 1.00 [ 1.00, 1.00] 0.916 NA NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.59*** [ 1.39, 1.85] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily persuasion utilized (partner’s view) 1.34*** [ 1.19, 1.54] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure experienced 0.95 [ 0.69, 1.29] 0.622 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pressure utilized (partner’s view) 1.49* [ 1.05, 2.29] 0.987 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing experienced 0.94 [ 0.71, 1.28] 0.665 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily pushing utilized (partner’s view) 1.28* [ 1.04, 1.64] 0.990 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Day 0.86 [ 0.64, 1.15] 0.843 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Own actionplan 9.48*** [ 7.77, 11.60] 1.000 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Partner actionplan 1.17 [ 0.96, 1.42] 0.941 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 1.40 [ 0.69, 2.84] 0.825 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean persuasion utilized (partner’s view) 1.33 [ 0.65, 2.67] 0.783 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure experienced 0.40* [ 0.18, 0.90] 0.986 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.51 [ 0.22, 1.15] 0.948 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing experienced 1.17 [ 0.41, 3.28] 0.616 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean pushing utilized (partner’s view) 2.07 [ 0.71, 5.84] 0.916 NA NA NA NA NA NA NA NA NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.30 [0.23, 0.41] NA 0.29 [0.23, 0.39] NA 0.60 [0.48, 0.78] NA 0.80 [0.46, 1.27] NA 1.15 [0.74, 1.76] NA
sd(Hurdle Intercept) 0.73 [0.56, 0.99] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily persuasion experienced) 0.11 [0.07, 0.17] NA 0.05 [0.02, 0.08] NA 0.04 [0.00, 0.10] NA 0.16 [0.01, 0.43] NA 0.22 [0.01, 0.51] NA
sd(Daily persuasion utilized (partner’s view)) 0.08 [0.04, 0.13] NA 0.05 [0.03, 0.09] NA 0.07 [0.01, 0.13] NA 0.21 [0.01, 0.52] NA 0.49 [0.13, 0.98] NA
sd(Daily pressure experienced) 0.06 [0.00, 0.22] NA 0.04 [0.00, 0.13] NA 0.07 [0.00, 0.24] NA 0.56 [0.11, 1.18] NA 1.08 [0.14, 2.43] NA
sd(Daily pressure utilized (partner’s view)) 0.06 [0.00, 0.18] NA 0.03 [0.00, 0.11] NA 0.08 [0.00, 0.26] NA 0.42 [0.02, 1.56] NA 0.83 [0.04, 2.75] NA
sd(Daily pushing experienced) 0.09 [0.01, 0.17] NA 0.07 [0.01, 0.15] NA 0.05 [0.00, 0.15] NA 0.21 [0.01, 0.51] NA 0.25 [0.02, 0.62] NA
sd(Daily pushing utilized (partner’s view)) 0.08 [0.01, 0.15] NA 0.03 [0.00, 0.09] NA 0.07 [0.00, 0.17] NA 0.15 [0.01, 0.62] NA 0.26 [0.01, 0.96] NA
sd(Hu Daily persuasion experienced) 0.21 [0.03, 0.40] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.19 [0.02, 0.37] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure experienced) 0.18 [0.01, 0.68] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.23 [0.01, 0.89] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing experienced) 0.57 [0.30, 0.98] NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.22 [0.02, 0.54] NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.68 [0.65, 0.70] NA 0.57 [0.56, 0.59] NA 0.96 [0.94, 0.98] NA NA NA NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()